Stochastic Gradient Langevin Dynamics
見やすいページを作ったのでこちらを見てくださいな
はじめに
目標: ある分布$ \pi(x)からのサンプリングを行いたい 1. Metropolis-Hastingsアルゴリズム (MH)
2. Hamiltonian Monte Carlo (HMC)
3. Langevin Dynamics (Metropolis-adjusted Langevin Algorithm)
4. Stochastic Gradient Langevin Dynamics (SGLD)
の順に見ていくと理解しやすい
Metropolis-Hastings
Metropolis-Hastingsについては既知のもとする
提案分布 $ q(z)を元に判定関数を用いて受容・棄却を行うMCMC 提案分布が対称, すなわち$ q(z|z_*) = q(z_*|z)ならば, 判定関数はただの$ \min\{1,\frac{\pi(z_*)}{\pi(z)}\}となる
これを単にメトロポリス法と呼ぶ
Hamiltonian Monte Carlo
Hamiltonian Monte Carlo (HMC)
$ H := H(q,p;t) =\mathcal{K}(p)+\mathcal{U}(q)
HMCのお気持ち
解析力学に則り, ある系における座標と運動量の時間に関する軌道を元に, サンプリングを行う
例えば, $ (z,p)があるとして, 一定時間が経ったときの終点$ (z_*,p_*)が得られればサンプル取れるよね
ハミルトニアン自体は系全体のエネルギーそのものなので, エネルギー保存則が成り立つから
簡単のため, 質量を$ 1としておく (下記のように$ \mathcal{K}が表せるため)
$ \mathcal{K}(p) = \frac{1}{2}p^\top \mathrm{I} p
サンプリングしたい分布$ \pi(z) \propto \tilde{\pi}(z)について, 補助変数$ pを導入して, $ \pi(z,p) := \pi(z) = \pi(z)\pi(p)と拡張する. また$ \pi(p) = \mathcal{N}(p|0,I)とし, $ \ln{\tilde{\pi}(z)} = - \mathcal{U}(z)とすると,
$ \pi(p) = \frac{1}{{(2\pi)}^{k/2}}\exp(-\frac{\mathbf{p}^{\top}\mathrm{I}\mathbf{p}}{2})
から,
$ \pi(z,p) = \pi(z)\pi(p) = \exp(\ln{\pi(z)} + \ln{\pi(p)})
$ \propto \exp(\ln{\tilde{\pi}(z)}-\frac{{p}^{\top}\mathrm{I}{p}}{2})
$ = \exp(\ln{\tilde{\pi}(z)}- \mathcal{K}(p))
$ = \exp(-\mathcal{K}(p)- \mathcal{U}(p))
$ = \exp(-H(z,p))
となる
ここで, 提案分布を$ \mathcal{N}(0,I)とすると, Metropolis-Hastings→メトロポリス法となるので,
1. $ p を$ \mathcal{N}(0,I)からサンプリング
2. リープフロッグ法で$ (z,p)から候補点$ (z_*,p_*)を取得 (ステップサイズ $ \epsilon・ステップ数$ L) とすれば, 判定関数$ A(z_*,p_*)は
$ A(z_*,p_*) = \min\{1,\frac{\pi(z_*,p_*)}{\pi(z,p)}\} = \exp(-H(z_*,p_*) + H(z,p))
となる. だが, そもそもハミルトニアンは系全体のエネルギーそのものなので時間変化せず, 常に$ H(z_*,p_*)と$ H(z,p)は一致するので, 判定関数$ A(z_*,p_*)は常に$ 1になる めでたしめでたし
(実際はリープフロッグの誤差の影響で棄却されることもある)
ここからはハイパラのお話
リープフロッグではステップサイズ $ \epsilon・ステップ数$ Lというパラメタが存在する
$ \epsilonが小さければ探索空間が小さく非効率だし, 大きすぎると棄却される
なので, 基本は
①$ L固定で$ \epsilon動かす
②$ \epsilon固定で $ L動かす
のどちらかだが, 後者は後者で計算コストが高すぎる
したがって, ここらへんは慎重に選択する必要がある
リープフロッグ
一旦ここでリープフロッグのおさらいをする
$ p_i (t+\frac{\epsilon}{2}) = p_i(t) - \frac{\epsilon}{2}\frac{d\mathcal{U}}{dz_i}
$ z_i(t + \epsilon) = z_i(t) + \epsilon p_i(t + \frac{\epsilon}{2})
一つの式にまとめると
$ z_i(t + \epsilon) = z_i(t) - \frac{\epsilon^2}{2}\frac{d\mathcal{U}}{dz_i} + \epsilon p_i(t)
Langevin Dynamics
Langevin Dynamics (LD)
LDは先程のHMCの特殊例である. 具体的には $ L = 1としたときのHMCである
LDにおいて, あるモデルのパラメタ$ \thetaに関して, 確率分布$ \pi(\theta)からのサンプリングを考えてみる
上のリープフロッグの式より
$ \Delta \theta := \theta_{t+1} - \theta_{t} = - \frac{\epsilon^2}{2} \nabla \mathcal{U} + \epsilon p_i(t)
$ = \frac{\epsilon^2}{2} \nabla \ln{{\pi}(\theta)} + \epsilon p_i(t)
$ (\therefore \ln{\tilde{\pi}(\theta)} = - \mathcal{U}(\theta))
ここで, データ$ (X,Y)について, $ \pi(\theta) \leftarrow \pi(\theta|Y,X)であるとき, ベイズより
$ \Delta \theta = \frac{\epsilon^2}{2} \nabla \ln{\pi(\theta|Y,X)} + \epsilon p_t
$ = \frac{\epsilon^2}{2}( \nabla \ln{\pi(Y|\theta,X)} + \nabla \ln{\pi(\theta)}) + \epsilon p_t
$ = \frac{\epsilon^2}{2}( \sum_n\nabla \ln{\pi(y_n|\theta,x_n)} + \nabla \ln{\pi(\theta)}) + \epsilon p_t
ただし, $ p_tは先程の説明通り $ \mathcal{N}(0,I)からサンプリング
なので, これはノイズを加えた勾配法と同じようになる
Stochastic Gradient Langevin Dynamics
Stochastic Gradient Langevin Dynamics (SGLD)
お待ちかねのSGLDだが, これは単純に上の式の第一項をサブサンプリングするだけ
なので更新幅は
$ \Delta \theta = \frac{\epsilon^2}{2}( \frac{N}{M} \sum_{n=1}^M \nabla \ln{\pi(y_n|\theta,x_n)} + \nabla \ln{\pi(\theta)}) + \epsilon p
これで, 棄却率がほぼ0 ∧ 計算量激減のMCMCが出来上がったことになる